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Abstract 

A mathematical model is developed to describe the sound emitted from an arbitrary 
point within a turbulent flow near solid boundaries. A unidirectional, transversely- 
sheared mean flow is assumed, and the cross-section of the cold jet is of arbitrary shape. 
The analysis begins with Lilley’s formulation of aerodynamic noise and, depending 
upon the specific model of turbulence used, leads via Fourier analysis to an expression 
for the spectral density of the intensity of the fax-field sound emitted from a unit 
volume of turbulence. The expressions require solution of a reduced Green’s function 
of Lilley’s equation as well as certain moving axis velocity correlations of the turbulence. 
Integration over the entire flow field is required in order to predict the sound emitted 
by the complete flow. Calculations are presented for sound emitted from a plug- 
flow jet exiting a semi-infinite flat duct. Polar plots of the far-field directivity show 
the dependence upon frequency and source position within the duct. Certain model 
problems are suggested to investigate the effect of duct termination, duct geometry, 
and mean flow shear upon the far-field sound. 


1 Introduction 

There is considerable interest today in developing a new generation of supersonic and sub- 
sonic civil transports that are relatively quiet despite their enhanced capabilities. To achieve 
this noise reduction, improved theoretical models are required to account for effects ignored 
in previous theories of aerodynamic noise. 

In his celebrated analysis of aerodynamically-generated sound, Lighthill [1, 2] rearranged 
the equations of motion for a compressible fluid to obtain the inhomogeneous wave equation, 
the basis of classical acoustics for stationary media, for density fluctuations p 

VV-c^ = -s( y ,<)- (>) 
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Here, V 2 = d 2 jdy 2 is the Laplacian, c 0 the mean speed of sound, t time, and s an acoustic 
source that Lighthill identified as a convecting quadrupole related to the Reynolds stress 
tensor. Thus, one of Lighthill’s major advances, called Lighthill’s acoustic analogy , is that 
the problem of sound generated aerodynamically can be replaced by an equivalent stationary- 
medium problem involving an acoustic source distribution that incorporates the dynamics of 
the actual flow. The theory met with considerable early success in predicting major aspects 
of the radiated sound field of cold jets, and even today the theory and its direct extensions 
find practical use. Both Goldstein [3] and Lilley [4] provide more recent perspectives on 
Lighthill’s theory. 

However, as demonstrated by Lush [5] , careful comparison of jet noise data with Lighthill’s 
theory revealed subtle but significant discrepancies, such as the frequency- dependence of the 
far-field directivity, which the theory was unable to explain. It was gradually recognized that 
though Lighthill’s analogy is exact, there are difficulties in estimating the acoustic source 
distribution when the radiation is refracted by the flow of the medium. 

Inclusion of refractive mean flow effects requires modification of the operator in the gov- 
erning equation. Lilley [6] derived a nonlinear moving-medium wave equation from the exact 
transport equations to include these effects. Lilley’s equation is considerably more compli- 
cated than that used by Lighthill, and this proves to be the equation’s major disadvantage. 
Mani [7, 8], in an effort to retain the simplicity of Lighthill’s equation while including convec- 
tion of the medium, proposed a two- region plug-flow model: (1) an inner region possessing 
a plug flow surrounded by (2) a stationary region. Goldstein, on the other hand, employed 
a linearized version of Lilley’s equation for unidirectional, transversely sheared mean flows, 
and showed that in both the low [9, 10] and high [11] frequency limits mean flow shear can 
significantly influence the directivity in free jets. In these investigations, both Mani and 
Goldstein modeled the sound source as a convecting multipole. 

Lighthill’s theory, based upon the free-space Green’s function, also failed to account for 
the diffractive and sound-generating effects of solid surfaces near the turbulent flow, including 
the effect of duct termination. Goldstein and Rosenbaum [12] accounted for these surface 
effects in a Lighthill-type formulation. They predicted the far-field directivity in terms of 
certain correlation scales of the turbulence for a localized, convecting source. Because this 
theory is based upon Lighthill’s formulation, it does not include the refractive effect of mean 
flow upon the transmission of sound. However, one of its major advantages is that it can 
associate particular source regions as being relatively noisy or quiet. 

On the other hand, Mani [13] and others, including Savkar [14], Munt [15, 16], and Cargill 
[17]) included mean flow effects in the context of the plug-flow type models for terminating 
ducts. Mani considered flat ducts while the other investigators focused upon cylindrical 
ones. Each of these investigations assumed a waveguide, propagating from minus infinity 
within the duct, served as the source of the sound that radiated to the fax-field. Due to 
their waveguide sources, these models fail to include the convective amplification expected 
of turbulence-generated sound. Furthermore, these waveguide models fail to provide insight 
into the relative contributions of certain regions of the flow to the overall directivity. 

To investigate the refractive effect of mean shear in the presence of solid surfaces, we here 
adopt the approach of Goldstein and Rosenbaum but begin with Lilley’s equation rather 
than Lighthill’s. We thus obtain a general expression for the spectral density of the far-field 
intensity due to a unit volume of turbulence within a jet. Then, due to the difficulty in 
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solving Lilley’s equation and to gain insight into this class of problems, we apply the theory 
to a simplified plug-flow jet problem, and thereby derive explicit expressions to calculate the 
directivity for a localized source. 

The structure of the paper is as follows: In the next section, we derive expressions for 
the spectral density of the far-field intensity for a localized convecting acoustic source that 
is located within a turbulent jet. The analysis includes the effect of duct termination and 
the model applies to unidirectional jets with transversely sheared mean flow of any cross- 
section. Beginning with a linearized version of Lilley’s equation for the acoustic pressure and 
following a Fourier transformation in time, the transformed pressure is expressed in terms of 
an integral involving a reduced Green’s function of Lilley’s equation and a transform of the 
localized source. Integration by parts is then required to transfer derivatives with respect to 
source coordinates from the source term to the Green’s function. The final expression for the 
transformed pressure is then given by an integration over the source region of the contracted 
product of a Reynolds stress tensor and a second-ranked tensor involving gradients of the 
reduced Green’s function and of the mean fluid velocity. An expression for the spectral 
density of the intensity is then formed and via a set of variable transformations re- expressed 
in terms of a fourth-order moving axis velocity correlation function. We subsequently focus 
upon the spectral density of the intensity due to a unit volume of turbulence, and assume 
that the fourth-order correlations can be adequately represented in terms of second-order 
correlations to derive a simplified form of the spectral density. This expression is further 
simplified by assuming the turbulence possesses axisymmetric symmetry. A final expression 
for the spectral density is obtained upon assuming the turbulence is isotropic. 

Section 3 is devoted to the applying the foregoing theory to a plug-flow jet that exits 
between two parallel plates. A formal solution to the Green’s function is obtained using 
the Wiener-Hopf method of analysis. The analysis is followed through to its completion, 
and sample plots of the far-field directivity are presented in section 4 . Two appendices 
supplement the Wiener- Hopf analysis. 


2 General Theory 

Below we consider a jet of arbitrary cross-sectional shape that exits a semi-infinite duct (see 
Fig. 1 ) and develop a general expression for the spectral density of the far-field intensity due 
to a convecting localized acoustic source located within the flow. This general expression for 
the spectral density is then simplified by application of additional assumptions relative to 
the nature of the turbulence. 

2.1 General Expression for the Spectral Density 

Let the total fluid velocity v(y,t) at any point y in the fluid be given by the sum U(y) + 
u(y, t). Here, U represents the mean velocity, u the fluctuating part with zero mean, and t 
time. To simplify the analysis, we assume a unidirectional, transversely sheared mean flow, 
at least in a local sense. Then if L denotes a unit vector in the direction of the jet and 
y = (i/i) 2/2, 2/3) a Cartesian coordinate system with origin at the end of the duct, the mean 
flow is of the functional form iiU(y2, 2/3)- According to Lilley’s equation as presented by 
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Goldstein [3], the pressure variable II(y,t) = (c„/cp) lnp/po (with p = p(y, t) the pressure) 
is then governed by the third-order wave equation 


Wt (2) 

where the acoustic source T(y, t) is given by the expression 

r(y - *» = (§ v • v • uu - 2 < w > ' & ■ “") ■ (3) 

Lilley’s equation as given above is derived from the exact equations of motion upon neglect 
of viscous dissipation and entropy fluctuations and upon linearization about the mean flow. 
We further assumed that the mean speed of sound Co is constant and that the turbulence 
is incompressible (V • u ~ 0). Appearing above is the heat capacity ratio Cp/c v , a reference 
pressure p 0 , a reference density p 0 , and the material derivative ^ + U^. To insure 

convergence of integrals that will be encountered in the subsequent development, we assume 
that uu is identically zero for \t\ > T, where T represents a large period of time that at the 
end of the analysis may be taken as infinity. 

We next suppose that pressure fluctuations p — po are everywhere much less than the 
ambient pressure p 0 . This permits the replacement of II in Lilley’s equation with (p— Po)/PoCo> 
where p<> is the mean density. In writing this, we have taken the reference pressure as 
Pq = p 0 CQ with Co = yJcpPo/cvPo, which is valid for an ideal gas whenever fluctuations in 
entropy may be neglected. 

Solid boundaries are assumed to be rigid so that n • v must vanish identically. Here, 
n represents a unit normal vector directed from the fluid phase into the solid. To derive a 
boundary condition upon the pressure, we consider the normal component of the momentum 
equation: 


n 


dx 

p(-^-+ v - Vv )+ y p 


= 0 


(4) 


The first term is identically zero due to the assumption the boundary is rigid. For surfaces 
located outside the region of turbulent flow, the acoustic approximation applies so that the 
second term, which involves the square of a perturbation quantity, is negligible. The second 
term can also be shown to be identically zero for surfaces adjacent to the turbulent flow if the 
radius of curvature R of the surface in the flow direction is infinite. Otherwise, the second 
term is given by the potentially large quantity ±pv 2 /R. To make progress, we assume that 
R is much greater than, say h, a characteristic distance between boundary surfaces (see Fig. 
1). Under these conditions, we may require the normal derivative of the pressure to vanish: 


dp 
dn 

In addition to the above requirements, the pressure must satisfy causality and the outgoing 
radiation condition. 



0 on solid boundaries 


(5) 
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To solve this system of equations, we introduce the Fourier transform P( y) of the pressure 
with respect to time: 

1 /*00 

P = (p- Po)' = 4= / (P - Po)^' <“ TO 

yj L'K J-oo 

Transformation of equation (2) with II replaced by ( p — po)/ po<^ gives 

(- \k + m4-) V ! P - (-it + - 2 (VM) • A VP = r‘(y), (7) 

dyi [ dyi \ dyi 

with k = u;/co the wave number and M(j/ 2> 2/3) = U /cq the local Mach number of the jet. 
(The Fourier transform of a quantity with respect to time is denoted by a superscript t.) An 
expression for the Fourier transform of the source can be obtained from equation (3): 

r‘(y) = -Po (-i k + M^-)V • V • (uu)‘ - 2 (VM) • ^-V(uu)‘. . (8) 

The pressure boundary condition on solid surfaces transforms simply to = 0. Finally, P 

is required to satisfy the radiation condition at infinity as well as causality. 

A formal solution to P can be expressed in terms of a Green’s function G as 

P(x) = - G(x|y)r‘(y) dy 

r o pi ~ 

= Po f G(x|y) (—ik + M -: — )V • V ■ (uu)‘ — 2(VM) • ^ — V(uu)‘ dy. (9) 
Jv [ oyi °y 1 

Here, V represents all the volume outside of solid boundaries, but effectively, because the 
source is localized, the integration need only be performed over that region for which P is 
nonzero. The Green’s function satisfies the same governing equations and conditions as P, 
except that a delta function replaces the original source term. Thus, G(x|y) is governed by 
the field equation 

(-i* + m 4-) W x G - (-ik + m4~?G\ - 2 (V X M) • 4~ v * G = -^ x - y)> ( 10 ) 

OX 1 OX\ J OX\ 

with n • V X G = 0 on solid boundaries. Finally, we require G to satisfy a radiation condition 
at infinity and causality. Here, V x = ^ represents the gradient operator with respect to the 
observer’s coordinates x, and y now represents the source coordinates. 

It is convenient to integrate by parts in equation (9) and thereby transfer all the deriva- 
tives with respect to y to the Green’s function. To assist in this process, we define three 
integrals, 

I a = / GV • V • uu dy (11a) 

Jv 

r \ 

I b = [ GM-?-V V-uudy (lib) 

Jv dyi 

I c = [ G(VM) • 7p-V - uu dy (11c) 

Jv dyi 
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so that P may be expressed as the sum 


P = p 0 (—ikl a + I b - 2 I e )* (12) 

Next, we display certain easily proved vector-dyadic identities that are valid for arbitrary 
scalar, vector, and symmetric dyadic functions of position S'(y), V(y), and D(y), respec- 
tively: 

SV • V • D = V • V • (SB) - 2V • (D • VS) + D : VVS 
V • (SW) = 25VV • V + (VS) • VV 

„ av 2 dv 1 s-v 2 dv,s „ 

Vl5 '^7 - d Vl ' V2 

V V D = V ■ (D • V) - D : VV. 

In the third identity, Vi and V 2 are two distinct vectors. 

Now consider integral I a in conjunction with equation (13a) with S set to G and D to 
uu. Integrate the identity over the volume V and apply Gauss’s divergence theorem to the 
first two integrals. The resulting surface integrals vanish at infinity because both the source 
and the Green’s function vanish in this limit. They also vanish along all solid boundaries 
due to the condition n • v = n • u = 0. This is easily seen for the second surface integral. 
That the first integral vanishes along solid boundaries can be seen via the second identity 
above with S set to G and V to u. It follows that I a is given by 

I a = f uu : VVG dy. (14a) 

Jv 

In a similar manner, it can be shown using the above identities that I\> and I c reduce to the 
volume integrals 

I h = - [ uu : Aw(GM) dy (14b) 

oy i 

I c = / uu : ^-V(GVM) dy 

Jv dy i 

= / uu:^- [V(GVM)] S dy (14c) 

Jv oyi 

In the latter expression, the superscript S denotes the symmetric portion of the dyadic; e.g. 
if a and b are vectors, the symmetric portion of the dyadic ab is (ab + ba)/2. 

Putting these results together yields the expression for the transform of the acoustic 
pressure 

P(x) = -po^(uu)‘ : G(x|y) dy (15) 

with the second- ranked tensor G given by 

G(x|y) = iitVVG + -J- {VV(GM) + 2 [V(GVM)] 5 } . (16) 


(13a) 

(13b) 

(13c) 

(13d) 
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In this expression, G = G(x|y) and M = M{y 2 , 2 / 3 )- 

Because only statistical information is available for the source, we introduce the spectral 
density of the intensity 7 w (x). In the far-field at location x, the spectral density is related 
to P via the expression 


Li*) = 


\P? 

2Tcop 0 ' 


(17) 


As stated previously, the source is presumed to have a duration of 2 T with T some large 
period of time. Introduction of the above expression for P then yields 


/„(x) = 2 ^ X X G(x|yO : (uV)'Kur : G*(x|y") dy' dy" 


(18) 


for the spectral density. The asterisk denotes the complex conjugate. 

We now use the fact that the Fourier transform of a convolution is proportional to the 
product of the Fourier transforms of the two components, and that, provided f(t) is real, 
the Fourier transform of f(—t) is the complex conjugate of the Fourier transform of fit). It 
is then easily shown that 

T roo 

(u'u') < (u"u") t ’ = — / u'u'u"u"e _,u;T dr (19) 

7 T J-00 


where the fourth-order, two-point, two-time velocity correlation tensor is given by 

1 T 

u'u'u"u" = f u(y', t)u{y', t)u(y", t + T)u(y",t + t) dt (20) 

11 J-T 

Here, single primes on velocity u denote position and time (y 7 , t) and double primes position 
and time (y", t + r). 

We now decompose this correlation tensor as 

u'u'u"u" = 72 4 (y', y" — y', r) + u'u' u"u", (21) 


with 72-4 the fourth-order velocity correlation tensor 

72. 4 (y', y"-y\ r) = u'u'u"u" — u'u' u /# u". 


(22) 


Note that due to the assumption that the process is stationary, the correlation tensor 
u'u' u"u" is independent of r. Thus, integration of this term over r in equation (19) will 
lead to the delta function <5(u>). Consequently, this term is possibly nonzero only at zero 
frequency, u> = 0. An argument can be made that the contribution of this term to the 
spectral density is identically zero (see, e.g., Goldstein and Rosenbaum [12]) even for u = 0, 
a frequency which is of no physical interest. In any case, we now ignore this contribution 
and write the spectral density in terms of 72. 4 (y' y" — y', r) as 

/ “ (x) = 2^ J1 L L 6(x|y,) ; Ui{y ’' y " ~ y ’’ t) : <iy ' dy " dr - < 23) 

The above expression can be simplified a great deal. We first note the assumption that the 
turbulence is stationary in time leads to a useful symmetry property of 72. 4 (y', y" - y', r). 
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Specifically, in equation (22) interchange primed and double primed quantities and replace 
r by — r. Due to stationarity, the resulting integral corresponding to equation (20) must be 
independent of a translation in time. Upon replacing the variable t with t + r, we obtain the 
result in Cartesian tensor notation 

ft 4 (y', y" - y'> T ) = Kiwi y', v, T ) = Kkiiji y", »? (1 \ ~ T ) 

= \Kiikliy, 1, r) + iw«a(y*, •»">, -r), (24) 

with 17= y" — y' = — separation vectors. This expression is then introduced into the 
above expression for the spectral density, and it is easily shown that the second term is the 
complex conjugate of the first, thereby leading to the expression 

/o,(x) = f°° [ f G(x|y') : ft 4 (y\ V, t) : G*(x|y' + *7 )e~ i “ T <fy' drj dr, (25) 

27TCo J-oqJvJv 

where denotes the real operator. 

It is helpful to introduce two further changes in the variables of integration. We introduce 
a new vector y, 

y - (yV . (2« 

which together with 77 replaces y' in the above integrals. With respect to the original vectors 
y 1 and y w , y represents an average of the transverse locations but the same location as y 1 
with respect to the component in the direction of flow. 

A further change in the variables of integration takes into account the convection of the 
source: replace ( tj , r) with (£, f) via the definitions 

£ = rj - \iU c t, (27a) 

f = r. (27b) 

After this change in variables, replace the dummy variable of integration f with r and so 
obtain the expression for the spectral density 

7 w (x) = I If G(x|y + (ii£i — £)/ 2 ) : 

2ttcq J-ooJvJv 

n A {y, e, r) : G*(x|y + (i^i + Z)J2 + ii U c T)e~^ r dy d '£ dr, (28) 

where 

^(y, r) = 72-4 (y + (ii^i - ()/2, ( + ii U c t, r). (29) 

is the moving axis correlation function. In the above expression for the spectral density, V 
represents the same physical volume (but in the new coordinates) as in the original variables. 

Now, according to the argument of Goldstein and Rosenbaum [12], the typical correlation 
length l of the turbulence in the moving frame is small compared to a typical dimension h 
of the duct boundary (see Fig. 1); i.e., / < h. On the other hand, the length scales that 
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appear within the Green’s function (see equation (10)) axe the wavelength k~ l , h, and h^\ 
where h ^ is a characteristic distance over which gradients of the mean velocity change. If 
the correlation length is also much less than the wavelength and h^\ the following argument 
cam be made to simplify the above expression for the spectral density of the intensity: First, 
by definition of the correlation length, if |£| > /, then 7t 4 (y, t) ~ 0. But, because l is 
much smaller than any length scale that appears within G, the Green’s function is nearly 
constant for 0 < |£| < l. Consequently, under these conditions, f in the arguments of the 
Green’s functions may be set equal to zero in the above expression for the spectral density. 
This permits us to obtain 

4,(x) = 2 ^* J J v ®( x |y) : [/ ^ 4 (y, r) d£ : G*(x|y + iif/ c r)e _la ' T dy dr , (30) 

where the omission of integral limits denotes the integration is to be performed over all space. 
We also expect the above step to be valid for thin shear layers. If a shear layer of thickness 
is much smaller than the correlation length ( h ^ <C l) and the source is located outside 
of and not too close to the shear layer, we expect the Green’s function in the far-field to be 
relatively independent of both of these length scales. In any case, upon calculation of the 
Green’s function and knowledge of the length /, the validity of this step can be evaluated. 

For design purposes, interest usually centers not in the prediction of the total far-field 
sound, but rather in ascertaining the regions of the flow that produce an inordinate amount 
of sound in the far-field. For this purpose, we let J w (x |y) denote the spectral density of the 
intensity at x due to a source at location y: 

/«(x) = / /„(x|y) dy. (31) 

J V 

From the above relations, we see that / w (x|y) is given by 

/„(x|y) = {6(x|y) : jT [/ Mv, «. r) d{] : G*(x|y + dr} . (32) 

To obtain the total spectral density in the far-field, the above equation would need to be 
integrated over all the noise-producing regions. 

Further simplification of the above equation is not possible without making additional 
assumptions concerning the nature of the turbulence. 


2.2 Reduction in Order of Turbulence Correlations 

Goldstein and Rosenbaum [12] suggest representing fourth-order velocity correlations in 
terms of second-order correlations, not only for r = 0 as Batchelor had done, but for r ^ 0. 
Batchelor’s original suggestion (and we now resort to Cartesian tensor notation), 


UjUjUfrUf = u'u' u" k u'( + u'ujj. u"u" + u'-u' k at r = 0, 


is based upon the assumption that the part of the joint probability of the velocity (with zero 
time delay) associated with the energy-bearing eddies is approximately normal. Goldstein 
and Rosenbaum extend this reasoning by saying that if the velocity correlations are separated 


9 



in both space and time, the central limit theorem suggests the joint probability distribution 
would be even closer to a normal distribution. In light of this, equations (22, 29) indicate 
the fourth-order correlation 7£ 4 (y, £, r) = y, £, t) can be written as 

n ijk i( y, t t ) = n ik n }l + n it n jk , (34) 

where the second-order velocity correlation is defined by 

£, r) = uyl (35) 

When this result is introduced into equation (32) and the symmetric property of G is 
considered, the spectral density may be written as 

Ju,(x|y) = fa f°° [ n ik n :U dr 1 . (36) 

7 TCq L J — oc J 

Here, we have suppressed the x and y dependence of Gjj(x|y) and introduced the abbreviated 
notation 

G% s G* (x|y + hU c r) (37) 

2.3 Axisymmetric Turbulence 

Turbulence measurements (see [18]) indicate that IZij is approximately an even function of r, 
at least for jet mixing regions. Furthermore, under the assumption of locally axisymmetric 
turbulence, Goldstein and Rosenbaum [19] contend the correlation tensor may be represented 

by 

7t,j(y, £, r) = A£ t £j + B6ij + C6n8ij + D(8u£j + £ Xj &)> (38) 

where A , B, and C are functions of y, r, and £ = |£| and are even functions of £ x . The 
coefficient D is also a function of y, r, and £ = |£|, but is odd with respect to £ x . Goldstein 
and Rosenbaum [12] show that this assumption leads to the spectral density expression 

•f w (x|y) = ( f e~' UT GaGjjS — (G 22 G 33 + G 33 G 22 HQ 23 — Q 22 ) + Gy jGj^Qij drl .(39) 

7TCo l J—oo 1 > 

following extensive algebraic manipulations. (This expression differs slightly from that ob- 
tained by Goldstein and Rosenbaum in that we are unable to apply a reciprocity relationship 
that is valid for a Green’s function of the Helmholtz equation.) Here, Q i: = Q ]t and 


Qu = 

/ (*?i - *? 2 ) * 

(40a) 

Q\2 = Ql3 = 

J (^12 + ^ 11 ^ 22 ) d£ 

(40b) 

Q22 = Q33 = 

J (^22 — ^ 12 ) d£ 

(40c) 

Q 23 = 

J (^22 — ^ 23 ) d£ 

(40d) 

5 = 

[ n\ 2 m 

(40e) 


These five integrals depend only upon four correlations: 72. xx , 1Z 22 , 7t X 2, and 7^23- 
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2.4 Isotropic Turbulence 

Considerable simplification of the spectral density expression can be obtain upon assuming 
the turbulence is locally isotropic in the moving reference frame. Then coefficients C and D 
in equation (38) are identically zero, and the expression for the correlation tensor reduces to 

My, (, r) = Miii + BSij, (41) 

where A and B are functions only of £ and r. In this case, define a correlation scale L(y, r), 

My, r) = A-, I Rl di , (42) 

M 1 

where Vu^ is the rms of the turbulent velocity at y. Then Goldstein and Rosenbaum [12] 
show that 

Qk = \(?) 2 L (43) 

for i, j = 1,2, 3, and 

S = ](5?)’l. (44) 

Clearly, it follows that the spectral density of the intensity may be written as 

«*> = L C e_1 " T + L H - (45) 


It should be recalled that (5,j is a constant with respect to r in this expression. 

Determination of the far-field directivity due to a source at y has thus been reduced 
to finding the Green’s solution G of equation (10) in the far-field and then forming G, as 
defined in equation (16). Then, together with appropriate knowledge of the turbulence, the 
spectral density of the intensity at x due to a source at y may be found from equations 
(36), (39), or (45) according to whether the turbulence in the moving frame has no special 
symmetry properties, is axisymmetric, or is isotropic, respectively. 


3 Application to a Plug-Flow Jet Exiting a Flat Duct 

As mentioned previously, the Lilley equation that governs the Green’s function that appears 
in the above theory has few known exact solutions, and so it is expected that certain ap- 
proximations will be required to make progress. In this section, we apply the above theory 
to a simple model problem. 
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3.1 Formulation 

One of the simplest relevant problems is that of a 2-D, plug-flow subsonic jet of thickness 
2 b exiting a flat duct into a resting medium (see Fig. 2). In this model, vortex sheets of 
effectively zero thickness separate the jet outside the duct from the two fluid layers at rest. 
The equation (10) that governs G can then be replaced by two simpler equations in which M, 
the local jet Mach number, is constant. Current boundary conditions on the Green’s function 
must be supplemented with conditions across the vortex sheets [3]; specifically, pressure and 
particle displacement must be continuous across the sheets. Furthermore, we must specify 
that the sheets are attached to the end of the duct walls as well as their departing slopes. 
Application of these latter conditions require that we revert to a kinematical description of 
the problem; i.e., we must recast the problem in terms of a velocity potential . 

The assumption of a plug flow within the jet also permits simplification of the expression 
for G in equation (16). Specifically, we obtain 

G = VV£G, (46) 


where C = ik + M^-. Application of this operator to the equation governing G together 
with the relation LS(x — y) = (i k — M j^)8(x — y) yields 

(V 2 x + k 2 )CG =0 for|x 2 |>6 

V 2 X LG + k 2 (l + =6(x- y) for |x 2 | < b, ( ' n 

following integration. 

We now regard LG as the Fourier transform of “pressure” in time, and seek to solve 
the corresponding problem for the velocity potential x 2 ). It is evident that <f> satisfies 

(compare with Mani [13]) 

0 V 2 x + k 2 )<j> =0 for |x 2 | > 6 

v 2 x <t> + k 2 (i + xalr) 2 ^ = 5 ( x ly) for 1**1 < b 

in the bulk fluids. Due to the relation between velocity potential and pressure within the 
jet, the above source term S must satisfy 

ikcop 0 {l + (iAf/k)^-]$(x|y) = S(x - y). (49) 

< 7Xi 

Along the rigid walls of the duct, <f> satisfies the zero gradient condition 

d<f>/dx 2 = 0 for xi < 0 and x 2 = ±6. (50) 


Continuity of pressure along the vortex sheets requires 

M-) = [1 + (51a) 

[1 + (iM / k)^]4>(xi, -b+) = (51b) 

for x-i > 0. (The plus or minus signs following the symbol b denote the addition or sub- 
traction, respectively, of an infinitesimal positive quantity such that evaluation is taken on 
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the appropriate side of the duct walls or vortex sheet.) Furthermore, the displacements 
tj(x \ , ±f>) of the vortex layers from the respective mean positions at x 2 — ±6 are related to 
the velocity potentials via 

?(*.. ±») = ( 52a ) 
[1 + (i M/k^Uxu ±b) = ji-2i!a^| w=±>iF (52b) 

for xi > 0. In addition to the above requirements, we require the sheets to be attached to 
the ends of the duct walls and to depart with zero slope (Kutta condition). Furthermore, 
we seek a solution that satisfies both causality and the radiation condition; the Helmholtz 
instability, however, precludes full satisfaction of the latter condition. 

3.2 Formal Far-Field Solution via the Wiener-Hopf Method 

The above equations for the velocity potential constitute a classical Wiener-Hopf problem 
due to certain boundary conditions prescribed on the semi-infinite duct (xi < 0) with other 
conditions prescribed on the semi-infinite vortex layer (xi > 0). 

Noble [20] describes in detail how such problems are solved. The basic idea is to assume 
the wavenumber k has a small imaginary component, which here is taken as positive. Fourier 
transformation of the governing field equations with respect to the spatial coordinate x\ 
leads to a set of differential equations that can be solved with constants of integration to be 
determined from proper application of the boundary conditions. Fourier transformation of 
the boundary conditions lead to functions in the complex a- space that are analytic in either 
the “plus” or “minus” regions; these are denoted by subscripts + and — , respectively. Here, 
a is the Fourier variable associated with X\. The plus and minus regions are assumed to 
possess a small region of overlap due to the small imaginary component of k. The goal of the 
Wiener-Hopf method is to eliminate unknowns in such a way so as to obtain functions that 
are analytic in the minus and plus regions on opposite sides of an equation. Then, because 
of the common domain of analyticity, the one expression must be the analytic continuation 
of the other. The two sides of the equation then represent an entire function in the complex 
a-plane. If the solution is unique, it can be shown that the entire function is identically 
zero for all a. Integration constants can then be determined, and the solution in the Fourier 
transform space is formally determined. Inversion then leads to the desired result. 

We define the Fourier transform of 4> as $ = $+ + with functions $+ and analytic 
in the plus and minus regions, respectively. These functions are related by the equation 

$±(a,x 2 ,fc|y) = f H(±x i)<f>(x u x 2 ; %)e lfco,;ri dx x . (53) 

J — CO 

Here, H(x ) represents the unit Heavyside step function. In the subsequent Wiener-Hopf 
analysis, we usually suppress the dependence of $ and similar functions upon k , and y, and 
further adopt the notation of Noble; i.e., we simply write $ + (a,x 2 ) = $+(<*) = $+(£ 2 ) = $+ 
according to whichever form is convenient provided there is no risk of confusion. 

Fourier transformation of equation (47) then leads to the set of equations 

4>" — k 2 ~f 2 <& = 0 for |x 2 | > b (54) 

ikcop 0 (l+Ma) ' di 
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with primes denoting derivatives with respect to x 2 . Scalars 7 and vj are defined by the 
relations 


7(0) = (a 2 -l) 1/2 

= (a-l) 1/2 (a+l) 1/2 ( 56 ) 

07(0) = [a 2 — (1 + Ma ) 2 )] 1 / 2 

= (l-M 2 f 2 (a-a L )^ 2 (a-a v f\ ( 57 ) 


where ai = — 1/(1 + M) and ecu = 1/(1 — M). We further define the branch cuts for 7 and 
for vj as il to ±00 and ±(1 T M) -1 to ±00, respectively, and require both functions to be 
positive for a —* 00 — iO. The common strip S and the plus and minus analytical regions 
in the complex a- plane, R( + ) and #(_), respectively, are displayed in Fig. 3 . If we assume 
that k = |fc|exp(i£) and no other singularities are present in the region, the strip may be 
described as the region between the two lines in the complex cx plane that pass through the 
points a = <*l and a = 1 with slope - tan 8 . If 8 < tt /2 as displayed in the figure, the plus 
region is to the right and above the line passing through a = 07. Similarly, the minus region 
is to the left and below the line passing through a = 1 . It can be demonstrated that ft(/;7) 
and $l(kzu) axe positive throughout S. If other singularities are present within the described 
strip, the thickness of the strip can be reduced as necessary to exclude these singularities as 
long as there is a common region of overlap. 

It is convenient to decompose $ as $ = + V, where \P(x 2 ) satisfies the homogeneous 

portion of equations (54, 55). On the other hand, the function V(x 2 ) represents a free space 
Green’s function within the jet in that it satisfies the source term in equation ( 55 ), but no 
specified conditions on the duct or vortex sheets. In particular, we take V as 

f 0 for x 2 > b 

V (x 2 ) = < iexp(ifcoyi— V2l) f or ~ < J, 

l 2k*wcoPo{l+M a ) 10r X2 - 1 ' 


Application of the radiation condition to 'P yields the result 


#(*2) = 4 


Ae-^ 

Be ~ kwx 2 + Ce kwX2 
De k ^ X2 


for x 2 > 6 
for |x 2 | < b , 
for x 2 < — b 


( 59 ) 


where integration constants A,B,C, and D must be determined from application of the 
boundary conditions along the duct and vortex sheets. (There should be no risk of confusing 
these four constants of integration with identical symbols that appear in connection with 
the turbulence correlation tensor in, e.g., equation ( 38 ).) 

We now adopt the method advocated by Noble on page 125 of his book to obtain the 
pair of Wiener- Hopf equations associated with the current problem. Let S (o) and £> (o) be 
defined as the sum and difference of $ evaluated just outside the duct and vortex layers: 


S (o) = $(+6+) + #(- b -) ( 60 a) 

D (o) = V(+b+) - #(-&-) ( 60 b) 

The notation S ^ and is similarly used to denote the sum and difference of V P evaluated 
just inside the duct and vortex layers, and S (V) and D {V) the sum and difference of V 
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evaluated exactly at 23 = b and x 2 = —b. Primes, e.g., denote the sum or difference of 
the function derivatives with respect to x 2 evaluated at x 2 = ±6. These sums and differences 
may be decomposed into plus and minus functions, e.g., + S'- . Lastly, we note 

that because V is given in (58), the quantities S± and D and the related primed quantities 
axe known in principle at this stage of the analysis. As is seen later, the decomposition of 
a known function into plus and minus functions can be accomplished via an integration in 
the complex a- plane. 

It is easy to demonstrate from these definitions and equation (59) the following identities 
for the ‘outer’ sums and differences: 


S {o) = (A + D)e~ kby (61a) 

£>(°) _ (A-D)e~ kln (61b) 

s (o)' = gW _ _ k ^ A _ D)e~ kh (61c) 

Z>(°)' = D ( ° Y = —k~f(A + D)e~ kln (61d) 


We obtained the extra relation in the latter two equations from the assumption the duct 

(oV ioY 

is rigid, i.e., equation (50) requires S_ = D_ = 0. From the above set of equations, it 


clearly follows that elimination of A and D leads to 

£i o) ' = -jfc 7 Z) ( °) (62a) 

D+ y = -Jfc 7 S (o) . (62b) 

A similar procedure for the inner fields yields 

5^ = 2(C + B) cosh kbw (63a) 

D (,) = 2(C - B) sinh kbw (63b) 

S (iy = 2kw(C -B) cosh kbw (63c) 

D w = 2kw{C + B) sinh kbw, (63d) 

find elimination of B and C gives 

S^' y = kwD ^ coth kbw (64a) 

D^' y = kwS ^ tanh kbw. (64b) 


The above conditions derive largely from the differential equations and application of 
the radiation condition. We now consider the remaining requirements of conditions imposed 
along the duct and vortex sheets. The rigid boundaries on the inside of the ducts require 


S& + = X>!?' + D? y = 0. 


(65) 


The effect of imposing continuity of particle displacement along the vortex sheets can be 
discovered by eliminating 77 in equations (52 a, b) and taking the half Fourier transform. 
Addition and subtraction of the results lead to the two relations 


(1 + Ma)S{° y = 5 (,)/ + S (V) ' 
(1 + Ma)D { ° y = D w + D {vy , 
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(66a) 

(66b) 



wherein we used equation (65). 

The final set of conditions is found by defining the half Fourier transform of the pressure 
jump across each of the duct elements: 


F( a,±b) = J 

CO Q 

{^(*1,±4±) - [1 + (iAf/ArW — ]^(x a , ±b^-)}e ikaXl dxi 

-00 OX \ 

(67) 

Integrate and take the 

sum and difference to obtain 



SP = S<°> - (1 + Ma) ($» + S( v) ) 

(68a) 


DP = D<°> - (1 + Ma) (D« + £>< v >) . 

(68b) 


Due to the continuity of pressure along the vortex sheets, equation (51), it is clear that 
F(a,±b) is analytic in Consequently, we have added the minus subscript to 5 (F) and 

DW above. This concludes the formulation stage of the Wiener-Hopf problem. Below we 
proceed to solve it. 

Equations (62a, 64a, 66a, and 68b) are seen to constitute one set of four equations 
involving the five unknowns D^°\ D^\ 5+^ , , and DP , and the knowns and 

S( v )' = -kwD^ v) . Equations (62b, 64b, 66b, and 68a) constitute another set involving the 

unknowns S^°\ S^\ Dp , D®' , and SP\ and the knowns and DW — -kwS^ v \ 
Each set of equations leads to a Wiener-Hopf equation, which when split supply the missing 
relation. 

Straightforward algebraic manipulations of the above relations lead to the following set 
of Wiener-Hopf equations: 


F(a)5+ )f = —kDP - k( 1 + Ma)D (V) [l + tanhfckn] 
Z(a)Dp = -kSP - k( 1 + Ma)S (v) [ 1 + coth kbw]. 


Kernels K(o) and Z(a) that appear above are defined by the relations 


K(a) = i + 
Z(a) = 1 + 


(1 + Ma) 2 tanh kbw 

w 

(1 + Ma) 2 coth kbw 

w 


(69a) 

(69b) 


(70a) 

(70b) 


To make progress, we need estimates of the various terms in the above equations as |a| — ♦ oo 
in the strip, which in turn correspond to certain physical properties as xi — ► 0. According 
to the full Kutta condition, displacement of the vortex sheet behaves as 77(11, ±6) ~ x/ 
as -4 0+. This requires that both S^' and D ( p behave as 0(a -5/2 ) as |a| -» 00 in 
the positive half plane. The requirement that pressure on the duct be finite leads to the 
conclusion that DP and SP are of order 0(a q ) as \a\ — ► 00 in the negative half plane with 
q<l. 

As noted by Munt [16] and others, the vortex layer in a linear model such as this has 
an inherent instability, which grows exponentially rapidly downstream. Associated with this 
instability are zeros of the kernels Y and Z. which are located respectively, say, at u ( P and 
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in the complex a plane. The exact location of these zeros depend upon k; for real 
positive k they are located in the region 3? a < 0 and 3a > 0. (Here, 3 is the imaginary 
operator.) A numerical investigation revealed these zeros were nearly stationary when 8 is 
varied over the range 0 < 8 < 7r/2; furthermore, they were always located in the second 
quadrant of the complex plane, which we henceforth assume to always be true. 

The next step in the Wiener-Hopf procedure is to split the kernel functions as Y (a) = 
K t -(a)K.(a) and Z(a) = Z+(a)Z-(a ) with the plus and minus functions regular and non- 
zero in the respective half planes. If 8 is nearly zero, one can deduce that Y± and Z± behave 
as 0(a 1/2 ) as |a| — ► oo in the respective half planes. However, if 8 = 7t/2, the respective 
instability zeros become part of Y+ and Z + . Then Y + and Z+ axe 0(a 3 ^ 2 ) and YL and Z- 
are 0(a -1 / 2 ) in the respective half planes as |a| — *• oo. 

To impose causality [16, 21], we set 8 = x/2 and assume the latter split of the kernels. 
Rearrangement of the two Wiener-Hopf equations then leads to the following: 


y + 4 0) ' -J+ = J- 


z + d ( + y -k+ = k.~ 


kD { P 

Y. 

kS i F) 

z. 


(71a) 

(71b) 


The new functions appearing above are defined via 


J(a|y) 


— Ar(l -|- Ma)D^(a |y)(l + tanh kbw ) 

Y^a) 

J+(a|y) + J_(a|y) 


and 


K(a\y) 


—k( 1 + Ma) i S’^ v ^(a|y)(l + coth kbm) 

Z^a) 

AT + (a|y) + ff_(a|y). 


(72a) 


(72b) 


Here we indicate for future reference the dependence upon source location y through the 
free-space Green’s function V , (x 2 |y)- 

The left-hand sides of the two equations appearing in (71) are analytic in the plus half 
plane, and from the previous estimates, together with the evident smallness of J± and K±, 
vanish in the limit as |a| — * oo. Similarly, the right-hand sides are analytic in the minus half 
plane and vanish as |a| — > oo. Due to the existence of an overlap region, the two functions 
must be analytic continuations of each other. Moreover, application of Liouville’s theorem 
leads to the conclusion that both entire functions are identically zero. Thus we obtain the 
results 


s+^Vly) = 

|y) = 


*M<*|y) 

y+ 

(73a) 

|y) 

(73b) 
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The inversion formula for x? > b then leads to the expression 

£G(x|y) = «*“. ( 74 ) 

47 r J r 7 \ Y+{a) Z+(a) ) 

with the contour T a straight line defined via T : a = ue 16 with —00 < u < 00 (see Fig. 3) 
and 8 = 7t/2. 

We obtain the causal result for k real and positive by analytic continuation of the above 
formula. However, as <5 is decreased from tt/ 2 to 0, the contour T crosses the instability zeros 
of Y+ and of Z+. A residue term must be added to the result of the above integration for 
each of the two instability zeros whenever the contour lies on the other side of the respective 
zeros; thus for <5 = 0, the above formula must be supplemented by two instability terms. If 
these instability terms are not added, the resulting expression will not satisfy the boundary 
condition on the duct. 

Our interest centers upon the far-field, and there the instability waves are exponentially 
large near the downsteam axis, particularly for polar angles 0 less than approximately 45 
(See Appendix A for an investigation of the zeros of the kernels, which are related to the 
questions of causality and instabilities) as measured from the jet axis. There are various 
physical mechanisms not represented in the present linear model that will moderate the waves 
in this region. To make progress, we ignore the instability waves and limit our attention to 
predictions of the spectral density of the intensity for values of 6 greater than either angle 
associated with the instabilities. For this limited region, equation (74) with 8 = 0 is valid. 

We can apply the method of stationary phase to evaluate the integral in equation (74) for 
large distances from the origin. Before doing so, however, it is useful to make the instabihty 
zeros explicit and to express the kernel factors in terms of factors that arise from integration 
along the real axis. First, consider once again 8 = 7 r /2 and make the instability zeros explicit 
by writing T(a) = (0 - uj, y) )(a - sj, y) )f (a) and similarly for Z. Here, sf, y) is a second zero 
that exists in the cut plane, and equals the complex conjugate of u{, r) when <5 = tt/ 2. This 
new kernel T(q) could be factored in a straightforward manner by integration over the 
vertical contour. (See Appendix B for a discussion of factorization of the kernels.) Now 
both zeros he in the minus half plane, and thus must belong to Y+. Hence we may write 
Y+(a) = (a - t4, y) )(a - s^ y) )T + (a) and K_(a) = VL (a). A similar factorization exists for Z. 

Next, suppose 8 = 0+. In this case, the contour F collapses onto the real axis, passing 
above the negative real axis and the branch point at a = — 1 , and below the branch point ^at 
a = 1 and the positive real axis. In the process of decreasing 8 from 7t/2 to 0, the zero at u 0 
remains almost stationary, but the zero 4 , y) that was originally located at u^ ] migrates to 
a — —2/M — iO; this zero thus remains part of the minus analytic plane. If we designate the 
factorization along this contour with tildes, i.e., Y (a) = F+(a)y_(o), it clearly follows that 
y + (a) = (a - 4 y) )F+(a) and H(o) = Yl(a)/(a - 4 V) )> wit h similar results for Z. 

To apply the method of stationary phase, we define polar coordinates (r, 0) (See Fig. 
4 .) via the relations x\ = r cos# and X 2 — b = rsin#, and initially limit 0 to the range 
0 < 0 < 7 T. Deform the contour by setting a = — cos(0 + ir) with —00 < r < 00 . The new 
contour is a branch of an hyperbola, and the deformation does not cross the branch points 
at cx. — il. Ignoring any pole or other possible branch point contributions, the method [20] 
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gives for kr — > oo 


£G(x|y) ~ poco 



«M«|y) 


+ 


K + (a\y) 


(a - u ( 0 y) )V+(«) (a - ui Z) )Z + (a) 


je iA:r 


(75) 


with a = - cos$. Explicit expressions for J+ and K + are derived in the next subsection from 
equation (72) with 7_(a) and Z-(a) replaced by Y-(a)/(a - u^) and Z-(a)/{a - Uq Z) ), 
respectively. 

We now consider possible residue and branch point contributions to CG in the far field. 
Aside from poles due to the instability zeros and the branch point at a = 1, the integrand 
is free of singularities in the upper plane. The possible contribution of instability zeros 
is discussed below. The branch point is not crossed in the deformation, and thus cannot 
contribute to the integral. It is shown in the next subsection that the only singularities 
of J + and K+ in the minus analytic region are simple poles at the zeros of cosh kbvj and 
sinh kbzu, respectively; K + also possess a simple pole at a^,. However, Y+ and Z + possess the 
identical respective poles, and hence these poles cancel. As in the upper plane, the branch 
point at a = —1 is not crossed in the deformation, and may thus be dropped from further 
consideration. The only remaining possible singularities of the integrand are due to zeros of 
Y+ and Z + in the lower plane. There are in fact an infinite number of these zeros, but each 
gives rise to an exponentially small residue term in the fax field. Prior investigators, e.g., 
Munt [16], have shown that the sum of these terms is negligible compared to the dominant 
contribution given above. 

Consider the possible contribution of the instability zeros to the above integral. In 
particular, we wish to precisely define the region of the instability wave. Set 

u^ y) = - cos(^ y) + i t£ V) ) (76a) 


and 



(76b) 


where v!q * and Uq^ are the locations of the instability zeros when <5 = 0. It follows that the 
instability wave may be neglected relative to the contribution in equation (75) provided 6 is 
greater than the larger of 6 and 6^ (which as stated previously are each approximately 
45 °). 

Equation (75) represents a formal solution to the far-field Green’s function CG that 
appears in the expression G = VV CG (equation (46)). Successful use of this solution to 
predict the spectral density of the intensity requires several more steps, which are described 
in the next subsection. 


3.3 Spectral Density for a Source Located Inside the Duct 

We here use the above formal result for CG to derive an explicit expression for the spectral 
density of the intensity of sound in the far-field due to isotropic turbulence located at various 
positions within the duct. The primary steps needed to complete this task are 

1. Factorize Y and Z 
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2. Decompose J and K 

3. Evaluate “Correlation Integral” 

By “correlation integral”, we mean the integration involving the correlation scale L(y, r) 
over the variable t that appears in equation (45). 

Factorization of the kernels is by far the most difficult of these steps. Provided a kernel, 
say, y {(*)■, is regular and nonzero within the strip, and it possesses the proper asymptotic 
behavior, the corresponding factor y+ is given [20] by the expression 


ln ^ + (°) 2ttJ- 


co— iO 
oo+iO 


In y{z ) dz 
z — a ’ 


(77) 


f or ft 0+. The path of integration here is above the negative real axis, below the point 
a, and below the positive real axis. All singularities along the real axis must be assigned 
to either the upper or lower planes. The path of integration is then modified to pass either 
below or above, respectively, the singularity as displayed in Fig. 5. To avoid distracting the 
reader at this point of the analysis, we consider properties of the kernels and their detailed 
factorization in two separate appendices. In Appendix A we investigate the zeros of the 
kernels, and in Appendix B we consider a practical method to factorize them. In any case, 
the kernel factors are most conveniently calculated numerically. We here assume the kernel 
factors are available, and move on to consider the decomposition of J and K . 

The first step in decomposing J and K is to demonstrate that these two functions are 
regular at the branch points of w, namely, ay and ay. To show this, form the sum and 
difference of U(5) and V(-b) to obtain S and respectively, and substitute the 

results into equation (72). Upon expressing the hyperbolic tangent and cotangent functions 
in terms of exponentials, we obtain 


J(a) 

K(a) 


(a — tto^) exp(ifco:y 1 ) sinh fctZ7y 2 
ikcopoY- (a) & cosh kbzu 

(a — Uq Z *) exp(ifco:j/i) cosh kwy 2 
lkcopoZ-(a) OTsinh kbzu' 


(78a) 

(78b) 


These expressions clearly show that the only branch point of J and K is that associated 
with the kernel factors Y_ and Z_, namely the point a = 1. 

Decomposition of these functions is accomplished via use of a formula similar to that 
used for factorizations. For J + we have 


1 roo— iO 

J + (a) = — / 

Z7Ti J — oc-f-iO 


J(z) dz 

■) 

z — a 


(79) 


with an equivalent formula for K+. For a source located within the duct (t/i < 0), the 
contour may be closed over the lower half plane, and the integral evaluated using Cauchy’s 
integral formula. Summing over the residues yields 


J + (a) 


— 1 


k>bcop 0 {l - M 2 ) 


£ 


(— l) n (an-( 1/ i ) - 4 y) ) ex P( ifca Ui/2)3h) sin 


(2w— l)ffy 2 
2b 


Y-{a ( Jy /2) ){a ( J (l/2) - a)(a { ) 


( 1 / 2 ) 


a) 
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and 


K + (a) = 


k 2 bcop 0 


— (cto } - u^ } )exp(iA:aS V) 


2Z-(a { 0 - ) )(ai~ ) - a) 


1 


1 - M* £ 


(-l) n (a4 ] - t4, Z) )exp(ifca4 ^yx) cos nny 2 /b 


Z_(ai ^)(ai * - a)(a4 1 - a) 


(-) 


(80b) 


Here, 


1, \ M 

a ” 2 ai + au ^~ i - M 2 


(81) 


represents the midpoint between the two branch points of tu. The set of points > represents 

the zeros of w sinh(/cfou) and of cosh (kbw) in the lower-half plane: 


_) = / d - 1 - M 2 ) - for p > kb/(*Vl=W) 

a — ~ ^ 2 tt 2 (1 — M 2 ) for v < kb/(iry/]~—~M 2 ) 


a , ' = 


(82) 


Here, v = 0, 1, 2, . . . for zeros of to sinh(fc£>o7) and v = |, |, |, . . . for zeros of cosh(fc6o7). 
Observe that Oq _ * = oli,. Values of CG in the far field can now be determined from knowledge 
of the kernel factors together with the above expressions for J f and K+. 

Due to the exponential dependence upon y\ in both J+ and K+, every term of G = 
WCG has this same dependence, and this fact along with the form of the turbulence 
correlation found by Chu [18] for L( y, t) permits an exact evaluation of the “correlation” 
integral. 

To make progress in the following, we make this dependence upon y\ explicit by writing 
VVJ+(a|y) and VV J R'+(a|y) as 


WJ + (a:|y) = ^B„_ (1/2) (a)exp(ifcQ[ 1 _ ) (1/2) j/ 1 ) 

n=l 

oo 

VVif + (a|y) = ^B n (a)exp(ifca[ l - ) yi)- 


(83a) 

(83b) 


n=0 


The B „(a) are independent of y : and are easily determined from the above equations. In 
particular, if we define the set of dyadics F„ as 


F„ = exp(— iArQ ( t 7 ) y 1 )fe 2 VVexp(ifcQ([ / -) y 1 )sin(i/7rt/2/fe) 
for v = n — ( 1 / 2 ) with n = 1,2, 3, . . ., and 

F u = exp(— ifca ( t 7 ) y 1 )& 2 VVexp(i/:a', -) yi)cos(i'7r3/2/&) 
for v = n = 0, 1,2,.. ., the B^a) may be expressed by the equations 

-iMf’ - u< z) )F„ 


B„(q) = 


2k 2 bf i copo(aQ ] - a)Z-(a { 0 } ) 


(84a) 


(84b) 


(85a) 
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for v = 0, 

-i(-l).-H./»»( air > - 

- kWcotoil - AP'iYJair' jia 1 ,-' - a)(a!r ) - a) 
for 1 / = n — (1/2), n = 1, 2, 3, . . and 

B .. = i(-l)-(a<-’ - 4 Z> )F, 

" kWcopotl - - a)(ai- } - a) 


(85b) 


(85c) 


for v = n = 1, 2, 3, To complete the definition of B„(a), we give explicit expressions 

for F„: 

F„ = — ijiiPfe 2 a[,~ )2 sin(i/7rj/ 2 /6)+ 

(iji 2 + i 2 i 1 )ifc&ct^ — ^i/7r cos(i/7rj/ 2 /6) — i 2 i 2 i/ 2 7r 2 sin(i/7ry 2 /6) (86a) 

for i/ = n — (1 /2) with n = 1, 2, 3, . . and 

F„ = -[ixii/: 2 fe 2 0'^'"^ cos(^7ry 2 /6)+ 

(iii 2 + i 2 ii)ifc&a[f^i'7r sin(j/7 ry 2 /f>) + i 2 i 2 i/ 2 7r 2 cos(i/7n/ 2 /f>)] (86b) 

for v = n = 0, 1, 2, It then follows from equation (75) that G may be written as 


f, ( k y /2 I tr / E”,iB._ ( „nexp(itoL(im!li) E",o B, exp(ifcqj, 

6 V (a - 4 y) )y + (a) (oc-u[ Z) )Z + (a) ) 

with a = — cos 9 for kr — ► oo. 

Let Q T be the so-called correlation integral. For convenience, we here define the complex 
conjugate of the correlation integral as (compare with equation (45)) 

/ OO _ 

e iwT G(x|y + i 1 I7 c r)L(y, r) dr , (88) 

-OO 

and following integration take the complex conjugate to get Q T . From the above expression 
for G, it is clear the typical scalar integral /„ is of the form 

I , „ = r e iwT ( 1+M ' 0, ^ )) L(y, r) dr. (89) 


Repeated integration by parts leads to 


L = 


1 


w 4 (i + Mecdr 5 ) 


,s: 


°° ,iu/r(l+M c Qi h d T 


d'i r 4 


(90) 


Chu [18] experimentally measured the two-point space-time correlations of both the tur- 
bulent velocities and the squares of these velocities in the mixing region of a 4-inch circular 
jet for low subsonic Mach numbers. The experiments showed that the fourth derivative of 
the relevant “self-noise” data could be fit well with analytical functions of the type 


d*L 

dr 4 


= Lq v) sech (u^t) cos(u>„,t), 


(91) 
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with Lq V \ u>f, uj, experimentally-determined constants. Unlike the conditions of Chu’s 
experiments, the present model assumes a planar jet within a two-dimensional space. Per- 
haps most significantly, turbulent fluctuations are here constrained to lie within the two- 
dimensional space. However, for lack of better information, we will use these results in the 
present theory. Introduction of Chu’s expression into equation (90) followed by a standard 
integration leads to 


L = 




2kfk*4(l + M c ai' ] ) 4 


sech 


7r[(l + M c a[ })k + k m ] 


+ 


sech 


2kj 

tt[( 1 + - fc«] 

2k f 


11 - 


(92) 


where kj = w//co and k, = u;,/co are wavenumbers corresponding to the respective frequency 
parameters. 

The above result leads directly to the expression for Q T : 


q t 



B n _ (1/2) / n _ (1/2) exp(ifcQ ( n I ) (1/2) yi) 

(o - u^)Y + (a) 

E” =0 B n / n exp(ifco4~)yi ) \ V 

(a - u { 0 Z) )Z + (a) )] 


(93) 


with a = — cos 6 for kr — > oo. 

Our explicit expression for the spectral density of the intensity in the far-field follows 
from equation (45) and the above results. Specifically, we find 


-L(x|y) ~ 



(94) 


in Cartesian tensor notation. The above expression is used in the next section to predict the 
directivity for selected values of frequency, Mach number, and source location. 


4 Discussion of Results 

We discuss various aspects of the theory just derived to give insight into the graphical results 
that follow. A small subsection then gives suggestions for future work. 

Each term in equation (80) corresponds to a different mode of wave propagation within 
the duct, and hence represents the contribution of that mode to the far-field Green’s function 
CG. According to equation (82), only certain will be real, and the number of real 
increases with frequency kb] these correspond to traveling waves within the duct. Note that 
the associated terms in equation (80) are not exponentially damped in the far-field expression 
for CG and similarly for I u . Also note that at sufficiently low frequencies, Oq - * is the only 
real this term corresponds to plane wave propagation within the duct. 

All other will have a negative imaginary part. The corresponding terms in equation 
(80) ultimately make exponentially small contributions of order exp(— o4 - ^|) to CG 
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and a similar small contribution to l w . Here, O' denotes the imaginary operator. Note that 
as a source is moved upstream within the duct, these factors decrease rapidly with distance, 
so that for a source located sufficiently far upstream, only terms consisting of real a[r\ be-, 
the traveling modes, will to contribute to the sum. On the other hand, for a source located 
near the end of the duct, several terms with complex w iU contribute to the far-field, 
but the magnitude of each successive term will decrease exponentially rapidly. These terms 
ran be associated with acoustic waves that cannot propagate within the duct due to their 
wavelength being too large. However, if the turbulent source is located sufficiently near the 
end of the duct, some of these waves may escape to the far-field. 

The case where the source is located at dead center (2/2 = 0) is of particular interest. For 
a source locate there, only symmetric ( v = n ) modes modes contribute to C,G in the far field. 
However, both symmetric and asymmetric ( 1 / = n - 1/2) modes still generally contribute to 
the far-field intensity because it depends upon W£G. 

We now use the preceding theory to predict typical directivities of the far-field intensity 
due to sound emitted from a unit volume of turbulence within the duct. Frequency and 
source position are the primary parameters varied in this study. Before proceeding, specific 
values are required for the turbulence correlation parameters kj and k» that appear within 
equation (91) for the turbulent source. According to Chu [18], the data are best fit with the 
choices 

bokj = 2.3 M e 

bftk* — 

for a cylindrical jet of radius b 0 and exit Mach number M e . We accept these relations even for 
our two-dimensional space, and consider 60/6 as a geometric ratio (the radius ratio ) with 
a value near unity that can be varied to some extent. Also appearing in our expression for 
the spectral density is the eddy convection Mach number M c . We here assume the relation 

M c = jM. 

Goldstein and Rosenbaum [12] suggested that it may be suitable to take M c as the negative 
of that given above to model sound emitted from the entrance of a jet as in an augmentor 
flap; for this purpose, they principally employed a radius ratio of 0.25. As might be expected, 
they suggest a radius ratio of unity to model the sound emitted from the rear of the flap. In 
any case, different physical geometries might suggest the use of different values of this ratio. 
For our purposes, we take M c as positive, and use 0.25 and unity for the radius ratio to test 
the sensitivity of results to this parameter. 

The remaining parameters requiring selection are the Mach number M, the frequency 
parameter ifci, and the source location (j/i, 2/2)- Chu’s data was collected for low subsonic 
Mach numbers, and so we set M = 0.3 in all our examples. We consider three different 
frequencies, low, low intermediate, and high: kb = 0.25tt (wavelength A = 86 ), 0.75tt (A = 
2.676), and kb = (A = 0.336). For the lowest frequency, c>4 _) is the only real For 

the next frequency, a ( ~1 is also real. For the highest frequency, the real correspond 
to v = 0, 1/2, 1, ..., 6. Based upon the full duct width, these frequencies correspond 
to Strouhal numbers (w/2tt) 26/I7 of 0.83, 2.5, and 20, respectively. Finally, we assume the 
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source is located at various vertical positions near the end of the duct. More precisely, we 
set yi = — 6 , and present results for a source located near the lower wall (y 2 = —0.8756), in 
the center between the two walls (t/ 2 = 0), and near the upper wall (y 2 = 0.8756). 

We plot a scaled spectral density of the intensity I w : 

128rfe / cg/ a ,(xly) 

7 IfTH)' 

We inserted a scaling parameter Q above so that the maximum value of 7 W equals 100, 
thereby aiding comparison of the various results that follow. Figs. 6, 7, and 8 each contain 
four polar plots of the scaled spectral density. The three figures correspond respectively 
to the low, low intermediate, and high frequencies mentioned above. The first three plots 
within each set correspond to a source located near the lower wall, in the center, and near 
the upper wall. For these plots, the radius ratio equals 0.25. For the final plot of each set, 
the radius ratio is set to unity for a source located at (—6, 0). Thus, comparison of second 
and fourth plots within each set reveals the sensitivity to the radius ratio. Finally, a dashed 
line drawn at an angle of 45° also appears within each graph. The area to the right of this 
line represents the region where the instability wave is exponentially large, but which we 
choose to ignore in the present solution. 

Fig. 6 displays directivities for low-frequency radiation as predicted by the current theory. 
All four directivity curves are oval in appearance with maximums occurring near 40°. Looking 
more closely at Figs. 6A, 6B, and 6C, we note that as the source moves from near the lower 
wall to near the upper, the angle of the maximum increases slightly from 36° to 39°, and 
the maximum value of the spectral density monotonically decreases only slightly: Q changes 
from 2.27 (10~ 3 ) to 2.07 (10 -3 ). From Figs. 6B and 6D, we note the shape of the directivity 
curve is nearly unchanged by increasing the radius ratio by a factor of four, but the value of 
Q more than doubled with this change from 2.18 (10 -3 ) to 5.15 (10 -3 ). 

Directivities for our low-intermediate frequency appear in Fig. 7. Despite an increase in 
the frequency by a factor of three, the gross features of these directivity curves are similar to 
those of the lower frequency. Bbwever, closer examination reveals that many of the trends 
discovered above no longer hold at this increased frequency. As the source location moves 
from near the lower plate to near the upper in Figs. 7A, 7B, and 7C, the angle at which 
the maximum appears decreases from 44° to near 40°. Moreover, the maximum values Q of 
the spectral density in these respective cases are 1.53 (10~ 3 ), 1.34 (10 -3 ), and 1.45 (10 -3 ); 
there is no longer a monotonic decrease in Q as the source moves from the near the lower 
plate to near the upper. For a source located near the upper plate, Fig. 7C, a lobe appears 
forming near 60° that is absent is all the previous results. Figs. 7B and 7D suggest that 
the radius ratio again has little effect upon the shape of the directivity curve. But as in the 
low-frequency case, Q is sensitive to the radius ratio; it more than triples from 1.34 (10 -3 ) 
to 4.53 (10 -3 ) when this ratio changes from 0.25 to 1. 

Our high-frequency directivity results are displayed in Fig. 8. As at the lower frequencies 
examined above, the major portion of the radiation in the far-field occurs at angles 0 near 
40°. For a source located near the lower plate (Fig. 8A), multiple peaks are clearly seen with 
the magnitude of each subsequent peak decreasing as 6 increases. The source located at the 
center (Fig. 8B) gives rise to a single major peak near 40°. A much smaller peak in the 


-L(x|y) 


(96) 
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directivity is observed at an angle near 75°. The source located near the upper plate (Fig. 
8C) gives rise to multiple peaks near 40°. These maximums again decrease in magnitude as 
6 increases, but the angular distance between them are considerably less than for the source 
located near the lower plate. A minor peak near 80° is observed. In contrast to the lower 
frequencies examined above, the angle of the primary maximum is here nearly constant (40°), 
independent of the value of y 2 - The values of Q for these curves are 2.55 (1CT 5 ), 1.22 (10 -4 ), 
and 2.81 (10 -5 ), respectively; the trend represented by these Q-values is different than in 
either of the two preceding cases. Figs. 8B and 8D show two nearly identical directivity 
curves. However, the value of Q for the case of unit radius ratio is 5.21 (10 -17 ), several orders 
of magnitude less than the value 1.22 (10 -4 ) found for a radius ratio of 0.25. The sensitivity 
to this ratio can be traced to the exponential decrease in the values of the hyperbolic secant 
functions that appear in equation (92), which define the integrals High sensitivity of Q 
to the radius ratio can be expected at high frequencies. 

In the limit of zero jet Mach number, the present theory leading to equation (45) for 
I w is equivalent to the expression found by Goldstein and Rosenbaum [12] for the spectral 
density of “self-noise” intensity in their Lighthill formulation. Goldstein and Rosenbaum 
went on to apply their theory to a two-dimensional plug flow jet as we have done, except 
they assumed the resulting waves propagated within a three-dimensional space in contrast 
to our assumption of a two-dimensional space. Due to this difference in dimensions, it is 
difficult to quantitatively compare the two theories. Figs. 13a and 13b of Goldstein and 
Rosenbaum give the directivity of self noise for a source located at (-b,0) for low and high 
frequencies, respectively. Compared to present results, the directivities of Goldstein and 
Rosenbaum vary only slightly with the angle of observation from the plane of the jet. 

In connection with experimental evidence of the downstream beaming of low frequency 
radiation [5], we note that Goldstein and Rosenbaum find finite directivities in the limit 
6 — ► 0. In contrast, the present theory indicates that for finite jet Mach numbers the 
directivity vanishes in this limit, independent of the frequency. This property of the present 
theory, which is no doubt related to the well-known refraction of waves as they cross a vortex 
sheet [22], can be deduced from equation (70) and the realization that both y+(a) and Z + (a) 
must possess the singular behavior (a + l) -1 ^ 2 as a = — cos & — * —1. Consequently, the 
inverse of these kernel factors and hence I w both vanish in this limit. However, as noted 
above, the present theory should reduce to the self-noise theory of Goldstein and Rosenbaum 
in the limit of zero jet Mach number. We should therefore expect that if the limit M —* 0 
is taken prior to the limit 0 — ► 0, a finite directivity will result. That this is indeed the case 
can be seen by noting that the pole of the Z(a ) kernel at a = ox = c>4 ^ moves to —1 in 
this limit. However, according to equation (85a), Bo will also have a simple pole at a = — 1. 
These two singularities will therefore cancel in Bo /Z + , thereby leading to the conclusion that 
the present theory predicts a finite directivity along the jet axis in the limit of zero jet Mach 
numbers. 

As noted above, the present model predicts a zone of relative silence near the jet axis for all 
frequencies, not just for high frequencies as is seen experimentally. Thus, the present model 
cannot predict the downstream beaming observed by Lush [5] for low-frequency radiation. 
This failure of our model problem has at least two plausible explanations. One possibility is 
simply that the region of downstream beaming is also the region of wave instability, and our 
solution, which excludes this instability wave, is not strictly valid in this region. According 
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to this view, the directivity of sound within the instability region cannot be known from 
the present model. A second possibility is that the present model problem, selected for 
its relative ease in solution, does not capture the important physics that gives rise to the 
downstream beaming, e.g., the refractive effect of mean-flow shear. Goldstein [9, 10] showed, 
e.g., that convecting quadrupoles of infinite lifetimes in a free jet give rise to downstream 
beaming in the low frequency limit; that model predicts that the directivity of pressure in 
this region is proportional to the gradient of the Mach number at the source location. In 
the present plug-flow model, there axe no gradients of the mean velocity except within the 
vortex layers, which are effectively of zero thickness. We also precluded sources from being 
located within such layers. Thus, this possible cause of downstream beaming was eliminated 
by selecting a plug-flow model. 

In this contribution, we have insisted upon satisfying the condition of causality, and this 
insistence lead to the instability wave and the related question of interpretation within the 
instability region. We note, however, that Dowling, et al. [23] have argued that causality 
need not be satisfied in problems of the present type. 

Future Work 

Our expression for the spectral density of the fax-field intensity due to a unit volume of 
turbulence, equation (45), can be applied to several different geometries and mean velocity 
profiles. In the present effort, we applied the expression to a two-dimensional duct from 
which a plug flow jet exits. Below we list possible extensions of the present work. 

1. The present study examined only one Mach number and a few frequencies and source 
positions. Predictions of the model should be extended to a broader range of parame- 
ters, and these results should then be compared in detail with prior models and relevant 
experimental data. 

2. Compare predictions of the present two-dimensional space problem with those calcu- 
lated using a non-causal Green’s function. 

3. The present plug- flow model is limited to a two-dimensional space. It is anticipated that 
extension of the model to sound radiation in a three-dimensional space while retaining 
the two-dimensional duct /jet geometry would give insight into the directivity due to 
jets passing through narrow slits. 

4. Adaptation of this plug-flow model to cylindrical ducts would be helpful. 

5. The inclusion of mean flow shear (in whatever geometry it is feasible) would represent 
a tremendous advance in current understanding. 

6. Results of these investigations should be implemented into code (e.g., the MGB code) 
for the design of actual jets. 
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5 Summary 

Starting from a linearized version of the Lilley equation, we derived an expression for the 
spectral density of the far-field intensity for sound emitted from a unit volume of turbulence 
located within a jet and near solid boundaries. The theory is valid for arbitrary unidirec- 
tional, transversely sheared mean flows. The theory was applied to a simple two-dimensional 
jet that exits a relatively long duct. Due to the complexity of the governing equations, a 
simple plug flow was assumed for the jet. The resulting equations were solved exactly in a 
formal sense using the Wiener- Hopf technique. Using existing turbulence correlation data, 
the spectral density was calculated and plotted for a wide range of frequencies and for sources 
located at various points just inside the duct. 
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Appendix A Zeros of Kernels 

Due to the need to impose causality, the location of zeros of Y and Z (defined in equation 
(70)) in the complex plane and their path as certain parameters change play an important role 
in the present analysis. Where possible, we here employ analytical techniques to investigate 
the zeros, and, where these fail, indicate suitable numerical procedures. Munt [16] gives an 
excellent discussion on this topic in connection with a related kernel for cylindrical jets. 

For the present kernels, perhaps the most useful information is obtained upon assuming 
k is pure imaginary, i.e., 8 = arg k = 7t/2. In this case, all poles of Y and Z must lie along 
the real axis. Therefore, we can employ the argument principle around a large contour in the 
upper or lower plane to determine whether any zeros are present. Furthermore, the easily 
proved general identity, 


r*(a, k) = -Y{a*, - k *), (A.l) 

is extremely useful. (For this discussion of the kernel zeros, it is useful to consider Y and 
Z as functions of both a and k.) If fc is pure imaginary, for any zero located in the upper 
complex plane, say, at a = Uq * , there must also exist a zero of Y in the lower complex plane 
at . Exactly the same reasoning applies to the kernel Z. 

We examined the change in arg Y and arg Z along a closed path just above the real axis, 
say from — R < 3fta < i? for large R, and then along the semi-circular contour following the 
path a = Re'^ with 0 < <j> < tt. One can deduce (See Munt [16] for details.) that the change 
in argument for both kernels is 27T, indicating the presence of one simple zero in the upper 
plane for each kernel and hence another one at the complex conjugate of this root. 

If the roots are located in a region in the complex plane where 23ft kbuj » 1, then to 
within an exponentially small order the hyperbolic tangent and cotangent functions may be 
replaced by unity. To this order of approximation, the two kernels become identical to each 
other and to the kernel found by Crighton and Leppington [21] for a single semi-infinite 
vortex sheet. Crighton and Leppington located this zero at 

u 0 = - cos( - + ir 0 ), (A. 2) 


where r 0 is the positive root of 

cosh to = [1 + V\ + M 2 ]/ \flM (A. 3) 

This zero is located in the second quadrant of the complex plane, and the actual zeros 
a = Uq') and a = Uq Z ' 1 of the two respective kernels are expected to be located nearby. Since 
there must be a location in the upper complex plane where 23ft kbw » 1, it follows that 
these are the sole zeros in the upper plane indicated by the argument procedure. 

Before leaving this case, we note that the above argument procedure did not capture all 
zeros of Y and Z in the cut plane, for some zeros are possibly located along the real axis 
between —1 and ai. The number of zeros in this region increase with the magnitude of kb. 
There are no other zeros along the real axis. 

The above zeros account for all zeros of both kernels in the cut plane when k is pure 
imaginary. For 6 = 7r/2, the contour of integration T, used for factorization of the kernels 
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(See Figure 3.) coincides with the imaginary axis. Thus, all the discovered zeros lie in the 
minus analytic plane, and are part of the plus kernel factors. It is useful to trace these zeros 
as 8 decreases to zero, and thereby determine whether they become part of the plus analytic 
region. The trajectory of these zeros can be determined by numerically solving the implicit 
equation 


da 

08 


§Y_ 
ds_ 

dY’ 

da 


(A.4) 


For all cases investigated, we found 

1. The zeros at a = Uq Y ^ and a = Uq Z ^ remained virtually stationary as 6 decreased to 
zero. Therefore, these zeros were crossed by F, and are the so-called instability zeros. 

2. The zeros initially located at s[> y) = u (Y) ’ and s ( 0 Z) = u ( 0 Z) for 8 = 7t/2 migrated to just 
below the real axis for 8 — * 0. In this limit one can deduce these zeros are located (to 
within an exponentially small error) at a = —2/M. These zeros remain fully within 
the minus analytic region. 

3. The zeros initially located between -1 and a^ appeared to migrate about the point 
a L so as to lie fully below the real axis when 6 = 0. These zeros thus also remained 
within the minus analytic region. 


We now suppose k is real and positive. Both Y and Z possess an infinite number of poles 
associated with the zeros of cosh kbw and sinh kbw, respectively, along the line 3? a = a; 
they also possess a finite number of poles along the real axis between ai and au- 

Now consider the zeros. For 23? kbw » 1, the hyperbolic tangent and cotangent func- 
tions may again be replaced by unity. Alternatively, for —23? kbw » 1, these functions 
may be replaced by —1. This latter situation occurs in the first and third quadrants relative 
to an origin at a = a. The same procedure (rationalization followed by factorization of a 
sixth-order polynomial) given by Crighton and Leppington [21] can be used to identify all 
roots satisfying either criteria. Only two of the six polynomial roots are approximate roots 
of our kernels, and these were identified above: the instability zero at u 0 and the root near 
a = -2/M - iO. 

The other possibility is that 2|3? kbw\ < 0(1). This region consists of a strip on either 
side of the line 3? a = d and the strip on either side of the real axis between ai and au- It 
appears that these zeros are interlaced with the poles that occur in these regions. 

To more closely investigate these zeros, we factor 7coshfc&ro from the denominator to 
obtain 


7 cosh kbw 


cosh kbw + 


(1 -|- Ma ) 2 7 sinh kbw 


w 


(A.5) 


The quantity in square brackets is suitable for use with the argument method to investigate 
zeros in the cut plane because it has the same zeros as Y , but none of the poles. Furthermore, 
the branch point singularities at a = ±1 are weaker than in Y, which aids implementation 
of the argument method near these points. 
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We n um erically calculated the change in argument of the quantity above in square brack- 
ets along paths surrounding extensive regions of the upper and lower complex planes to 
determine the precise number of enclosed zeros. The results confirm that the only zeros in 
the cut plane consist of the instability zero and the zeros associated with poles in the strips 
mentioned above. The zero near a = -2/M is not captured because it lies on the branch 
cut. In addition, a finite number of zeros lie along the real axis between 1 and au', these 
are interspersed with poles of Y that he in this region. Numerical calculations indicate that 
these zeros move into the other Riemann sheet as 6 increases from zero. However, we have 
not been able to prove in general that this is the case. 

These results suggest that each kernel has only one zero, namely, for the Y kernel 
and u { 0 Z) for the Z kernel, that represents an instability zero. The contour of integration T 
does not appear to cross any other zero when S varies in the range 0 < 8 < x/2. However, 
since we have not been able to analytically prove this result, an investigation of the zeros of 
the kernels should be conducted for each new value of kb and M. 


Appendix B Numerical Factorization of Kernels 


Factorization of kernels tends to be the most difficult aspect of the Wiener-Hopf analysis, 
particularly if the goal is to obtain analytical expressions. However, the present analysis only 
requires knowledge of the kernel values at particular points in the complex plane, and for this 
purpose numerical evaluation of the factors is quite acceptable. In this appendix, we discuss 
a straightforward and practical approach for the numerical factorization of these kernels. 
Specific details are given for the factorization of the kernel Y. Factorization of the kernel Z 
can be treated in a similar manner. We first suppose a lies within the range — 1 < a < 1, 
and then later consider the case of complex a. 

We begin with equation (77) for the definition of the plus analytic factor TV ( ^ ) f° r the 
generic kernel y = 3V-T-- In using this expression, we assume no zeros exist in the strip 
and that y — i ► 1 as a — ► oo in the strip. Neither Y nor Z, defined in equation (70), satisfy 
this latter condition, but it is easy to multipy by suitable functions that have known factors 
and thereby correct this defect. For example, if we define y by 


Wl - M 2 
M 2 7 

/ 1 (1 + Mq ) 2 tanh \/l — M 2 

i 7 w J M 2 7 ’ 


(B.6) 


the condition is satisfied. The function 7 can be factored as 7+7_ with 7+ (a) = \A* + 1 and 
7_(a) = \A* — 1- Then suitable factors of Y = Y+Y~ are V+(o) = M 2 y+(a)y + /\/l — M 2 
and Y-(a) = 7 _(a)T_. (Recall that we use a tilde to denote the factorization Y = Y+Y- 
along the real axis, and similarly for Z.) 

Equation (77) derives (See Noble [20]) from application of Cauchy’s integral theorem 
to a thin rectangular domain that lies within the analytic strip S. The path of integration 
represents one of the long sides of the rectangle. Because our interest assumes k is real and 
positive, the limit arg k — ► 0 must be taken prior to numerical evaluation. This limiting 
process causes the kernel y, which by definition is analytic within the strip (of finite width 
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when argJfc > 0), to appear to have singularities along the path of integration. We have 
to deal appropriately with these singularities. In particular, it is important to know the 
locations of the singularities, their type, and which side of the singularities the path of 
integration follows. 

Over the path of integration, different branches of the logarithm function are normally 
encountered. It is necessary to keep track of these branches during the integration process. 

We note certain easily-determined properties of y along the path of integration (See 
Figure 5) for k real and positive. First, y is real, positive, and free of singularities for 
a < — 1 and for a > ay. Simple poles due to zeros of j 2 are present at a = ±1. Also, 
depending upon the values of kb and of M, cosh kbtu may vanish along the real axis between 
ax, and ay, thereby giving rise to 2 N p simple poles of y in this region. The location of these 
poles to the left of a are given by equation (82) for v = n — (1/2), n = 1 , 2, 3, . . . , N p 
with v < kb/xy/i - M 2 . The location of poles to the right, say a[, +) , can be obtained from 
c*(/) = 2d — ajf). There also may be zeros located in the range 1 < a < ay. The locations 
of these zeros can generally be determined using a secant method; they alternate in order 
with poles that are present. We designate these roots as Uj with j = 1, 2, 3, . . . , N z . It is 
also easily deduced that the imaginary part of the kernel is identically zero in this range. 

We know the path of integration lies above the singularities at —1 and o^~\ and below 
those at a, 1 and a[, +) . To determine whether the N z zeros lie in the upper or lower planes, 
we use equation (A.4) to determine the signs of the derivatives da/d6 for <5 = 0. In all cases 
studied thus far for both kernels Y and Z, this sign has been positive, indicating each zero 
belongs to the upper analytical region. Thus, we here assume the path of integration goes 
below these zeros. 

As observed above, y is real, positive, and regular outside the range — 1 < a < ay, and 
hence the (numerical) integrations in these regions is straightforward. However, within the 
range — 1 < a < ay, a suitable representation for In}* is required, and to aid us in this 
task it is useful to define a related function, say, ip s . This new function is defined such that 
arg^e - '^ 5 ) is continuous at the singular points of y. (Compare with Munt [16].) Here, the 
change in arg >’ as the point a moves from the left to the right of a singularity is always ±7r 
since only simple zeros and poles are present. Whether the argument increases or decreases 
depends upon whether the singularity is a zero or a pole and whether the path goes above 
or below. Using our above knowledge of the singularities, we define t/>s via the equation 

^s(<*) = 7T {[H(a - 1) - H(a + 1)]+ 

£[tf(a - ai-a/2)) -H{a- afc^^)] - £ H{a - Uj )} 

n=l j = 1 


or, more conveniently, 


N S 

= * (B.7) 

3 - 1 

Here H(x ) is the unit step function, N s = 2 N p + N z + 2 the total number of singularities, 
qj the coefficient (either ±1) that corresponds to the jth singularity, and A j the location of 
the singularity. 


34 



We next define a function, say that takes into account changes in the branch of the 
logarithm function. This function has a form similar to that of ips : 

N b 

*,(«) = H(a — Xj), (B.8) 

i=i 


Here Nb is the number of such jumps and Xj the location of the j th jump. The coefficient 
qj is ±1 and corresponds to whether the argument of y is increasing (+) or decreasing (-) 
at Xj. The locations of branch jumps a = Xj are found by numerical solution of 

argMaje^-^-^U^ = 0, (B.9) 

for j = 1, 2, 3, . . . , Nb- Here, rp a , which we may restrict to the range 0 < xp a < ir, is a 
constant angle that may be used to adjust the location of branch jumps away from the Xj 
singularities. Setting this quantity to 7 r has been useful in the current investigation, but 
other values may be required for other Mach numbers and frequencies. 

For a in the range — 1 < a < a v , we now write ln^V as 

ln^(a) = {lnp^lcvje 1 ^-^^))] + i^H) + i[V>s(<*) - 4>a], (B.10) 


where ln p z denotes the principal branch of In z with — t < arg z < t. The imaginary part 
of the quantity in curly brackets is continuous in the given range of a, and the real part has 
only logarithmic singularities at the zeros and poles represented by the Xj. 

To proceed with the determination of y+(a), we evaluate the half-residue term in equation 
(77) to obtain the expression 




(B.ll) 


where V denotes the principal value. Noting a key result of analytical integration, 


1 v j au i ^s(^) dz 
27ri 7-i z — a 


1 (N S \ : n s 

2 IS 1i I ln i>( a u - - 2 H 1) ln » A J - “I. 


(B.12) 


we find that 34(a) can be efficiently calculated from the following equation for a in the 
range — 1 < a < 1: 


In 34(a) = ^{hi p [;y(a)eW»-^] + i[M<*) + M<*) ~ 

1 y- 1 ln p y(z) dz _\_p r a u {ln p [y(ar)tf^«~^(*M] + iij; B (z)} dz + 

27ri J-oo z — a 27ri J-\ z — a 

\ (| «) M<* - a) - i| * to, |A, - «l ■ - £ 1% ) 

1 r°° ln p y(z) dz 

2 TTi Joty Z Of 
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The above equation is not valid for a = Xj, but any needed values of y + (Aj) can be obtained 
by evaluating the above expression for nearby values of a and interpolating the results. 

The principal value integral in the above expression is evaluated numerically as follows: 
Constants Xj are placed in consecutive order to form subintervals, and integrations are 
performed over these individual ranges. A numerical method (e.g., the Gaussian-Kronrod 
adaptive method) that does not use endpoint evaluations is required. Only the integral over 
the range containing the point a needs to be calculated as a principal value. 

The numerical integrations in the above formula can be time consuming, and so it is 
convenient to represent the functions Y+ and Y- using splines, polynomials, or rational 
functions for real a within their respective domains, a > a and a < a, of analyticity. 
Outside these respective domains but within the range — 1 < a < 1, the factors can be 
obtained from the relation Y = Y + Y- . 

For our directivity calculation, we also need K_(a) and Z_ (a) for a = in the lower 
complex plane. These can be determined from equations of the type 


h y_ (a) » 

2 xi J- 00+io z — a 


(B.14) 


where the path of integration is the same as above except it is not deformed around a; i.e., 
there is no need to calculate the residue term or the principal value. Straightforward analysis 
shows that for a in the lower complex plane, T_(o:) may be calculated from 


ln3?_(a) = 

1 /•-» b p y(z) dz 1_ r a v {b p [y(z)efofr«~* g t*M] + iV> B (z)} dz _ 

27ri J — oo z — a 2ni J-i z — a 

1 \b 

2 II hl la p{ a v - a) ~ - a)] + ^[ln p (ay - a) - ln p (-l - a)] - 

J_ /“ (B.15) 

2tt\ Jau z — a 

As in the case for real a, the integral over the range — 1 < z < au is to be subdivided 
into several integrals so that no pole or zero of the kernel is crossed during the course of 
evaluation. 
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Fig. 1 : Side view of jet with arbitrary 
mean velocity U. Convecting source 
displayed. 



Fig. 2: Two-dimensional duct and plug-flow 
jet for Wiener-Hopf problem. Vortex sheet 
displayed. 
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Im(a) 



Fig. 3: The common strip S in 
the complex a-plane with plus and 
minus analytic regions indicated. 
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Fig. 4: Definition sketch of polar coordi- 
nates (r, 0). Also shows coordinates for 
moving eddy source at y = (y l5 y 2 ). 

(+) 



Fig. 5: Sketch of complex a-plane for 
factorization of Y showing contour of 
integration around singularities. 
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Fig. 6: Polar plot of scaled spectral density of the far-field intensity as a 
function of polar angle 0 for frequency kb= 0.25 n, exit Mach number M= 
0.3, convection Mach number Mc=0.15, yi/b=-1.0, and (A) y 2 /b=-0.875 and 
radius ratio bo/b=0.25; (B) y 2 /b=0 and bo/b=0.25; (C) y 2 /b=0.875 and 
bo/b=0.25; and (D) y 2 /b=0 and bo/b=l . Dashed line shows region of 
instability wave. Q-values for the respective plots: 2.27 (10' 3 ), 2.18 (10' 3 ), 
2.07 (10‘ 3 ), and 5.15 (10' 3 ). 
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Fig. 6: Continued 
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Fig. 7: Polar plot of scaled spectral density of the far-field intensity as a 
function of polar angle 0 for frequency kb= 0.75 n, exit Mach number M= 
0.3, convection Mach number Mc=0.15, yi/b —1.0, and (A) y 2 /b— 0.875 and 
radius ratio bo/b=0.25; (B) y 2 /b=0 and bo/b=0.25; (C) y 2 /b= : 0.875 and 
bo/b=0.25; and (D) y 2 /b=0 and bo/b=l . Dashed line shows region of instability 
wave. Q-values for the respective plots: 1.53 (10' 3 ), 1.34 (10‘ 3 ), 1.45 (10 -3 ), 
and 4.53 (10’ 3 ). 
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Fig. 8: Polar plot of scaled spectral density of the far-field intensity as a 
function of polar angle 0 for frequency kb= 6 tc, exit Mach number M= 0.3, 
convection Mach number Mc=0. 15, yi/b =-l .0, and (A) y2/b=-0.87 5 and 
radius ratio bo/b=0.25; (B) y2/b=0 and bo/b=0.25; (C) y2/b=0.875 and 
bo/b=0.25; and (D) y 2 /b=0 and bo/b=l. Dashed line shows region of instability 
wave. Q-values for the respective plots: 2.55 (10‘ 5 ), 1.22 (10 4 ), 2.81 (10' 5 ), 
and 5.21 (10 17 ). 
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